AF wavefunction differential imaging dev¶
Exploring differential imaging, aligned vs. unaligned samples. Currently computed via effective AF wavefunctions, expanded in spherical harmonics (see Formalism below), but may not be the best way (should be able to do via matrix elements directly?).
For final phase results, see Comparison with expt section.
22/09/23
Tidy version for distro - final results only.
07/09/23
Mainly working from previous AF notebooks, see http://jake/jupyter/user/paul/doc/tree/projects-share/carlos_tross_HHG_phases_2023/notebooks/AF_wavefns_method_dev_010221_class.ipynb
_basic-demoversion, for N2 with test alignment only. (Working with 3d-AFPAD-dev branch, requirespadPlot()fix as per this commit, testing on Jake inepsdev-shared-100122kernel.)
NOTE: currently using ePSproc base class, should update to PEMtk base.
Setup¶
# Imports
import numpy as np
import pandas as pd
import xarray as xr
# Special functions
# from scipy.special import sph_harm
import spherical_functions as sf
import quaternion
# Performance & benchmarking libraries
# from joblib import Memory
# import xyzpy as xyz
import numba as nb
# Timings with ttictoc or time
# https://github.com/hector-sab/ttictoc
# from ttictoc import TicToc
import time
# Package fns.
# For module testing, include path to module here
import sys
import os
# modPath = r'D:\code\github\ePSproc' # Win test machine
# modPath = r'/home/femtolab/github/ePSproc/' # Linux test machine
# sys.path.append(modPath)
import epsproc as ep
# TODO: tidy this up!
from epsproc.util import matEleSelector
from epsproc.geomFunc import geomCalc
from epsproc.classes.base import ePSbase # Data class
from epsproc.classes.multiJob import ePSmultiJob
# Plotters
import matplotlib.pyplot as plt
from epsproc.plot import hvPlotters
hvPlotters.setPlotters()
# Skip warnings for old XR versions
import warnings
warnings.filterwarnings('ignore', category=FutureWarning) #message='model will result in indexing errors*')
Load data¶
Testing for $N_2$.
# Load data from modPath\data
from pathlib import Path
modPath=Path(ep.__path__[0]).parent
dataPath = os.path.join(modPath, 'data', 'photoionization')
dataFile = os.path.join(dataPath, 'n2_3sg_0.1-50.1eV_A2.inp.out') # Set for sample N2 data for testing
data = ePSbase(fileIn = dataFile, verbose = 1)
# data = ePSbase(dataPath, verbose = 1)
data.scanFiles()
# N2O data
# dataPath = r'D:\VMs\Share\ePSshare\N2O\N20_wf' # Test dir on Stimpy (Win machine)
# data = ePSmultiJob(dataPath, verbose = 0)
# keys = [0,1,2] # Set for 1s datasets only (ugly!)
# data.scanFiles(keys = keys)
# # Scan data file
# dataSet = ep.readMatEle(fileIn = dataFile)
# dataXS = ep.readMatEle(fileIn = dataFile, recordType = 'CrossSection') # XS info currently not set in NO2 sample file.
data.jobsSummary()
Matrix elements¶
# key = 'orb1' # Set for O1s
# key = 'orb2' # Set for N1s (central)
# key = 'orb3' # Set for N1s (terminal)
data.calcOpts['thres'] = 1e-4 # TODO: propagate global threshold setting! Used for calcs but not for plots at the moment.
# Plot matrix elements
# data.lmPlot(keys = key)
data.lmPlot()
Default case, unaligned¶
In this case the axis distribution is isotropic, and there is only a single term with $A_{Q,S}^{K}=A_{0,0}^{0}=1$.
We'll calculate all the terms in the expression for $^{AF}T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\hat{k}_{L})$, but skip the summations and multiplication by $Y_{l\Lambda}^{*}(\hat{k}_{M})$, thus defining effective LF/AF matrix elements.
# Default alignment terms,
AKQS = ep.setADMs()
# ep.lmPlot(AKQS, xDim='t')
AKQS
# Run default af numeric code
# This will default to Type='L', it=1 and an isotropic axis distribution.
data.afpadNumeric()
# Set as ref data
key = 'orb5'
for item in ['TXaf', 'DeltaKQS']:
data.data[key][item + 'Unaligned']= data.data[key][item]
With alignment¶
Here we can look at the dependence of the channel coupling parameters $^{AF}\Delta_{l,m}(K,Q,S)$, as well as the final results for a given axis distribution...
Channel couplings for various $K$¶
In this case, assume a $z$-polarized ionizing field, and a cylindrically symmetric distribution. To look at the channel couplings, set all terms $A_{Q,S}^{K}=A_{0,0}^{K}=1$ for even-$K$ up to $K_{max}=8$.
# # Set ADMs - TEST ONLY
# KQSlist = [[K,0,0,1] for K in range(0,9,2)]
# AKQS = ep.setADMs(ADMs = KQSlist) # TODO: routine to automate this type of setting!
# # daPlot, daPlotpd, legendList, gFig = ep.lmPlot(AKQS, xDim = 't', pType = 'r', squeeze = False) # Note squeeze = False required for 1D case (should add this to code!)
# # daPlotpd
# AKQS #.coords
# SET ADMs from file, see https://pemtk.readthedocs.io/en/latest/fitting/PEMtk_fitting_basic_demo_030621-full_010922.html
# Load time-dependent ADMs for N2 case
# Adapted from ePSproc_AFBLM_testing_010519_300719.m
from scipy.io import loadmat
ADMdataFile = os.path.join(modPath, 'data', 'alignment', 'N2_ADM_VM_290816.mat')
ADMs = loadmat(ADMdataFile)
# Set tOffset for calcs, 3.76ps!!!
# This is because this is 2-pulse case, and will set t=0 to 2nd pulse (and matches defn. in N2 experimental paper)
tOffset = -3.76
ADMs['time'] = ADMs['time'] + tOffset
AKQS2pulse = ep.setADMs(ADMs = ADMs['ADM'], t=ADMs['time'].squeeze(), KQSLabels = ADMs['ADMlist'], addS = True)
# Subselect
# AKQS2pulseSub = AKQS2pulse.sel(t=slice(4, 5, 4))
AKQS2pulseSub = AKQS2pulse.sel(t=slice(4, 6.5, 4))
# AKQS2pulseSub
# PEMtk class has this directly:
# data.setADMs(ADMs = ADMs['ADM'], t=ADMs['time'].squeeze(), KQSLabels = ADMs['ADMlist'], addS = True)
# data.data['ADM']['ADM']
# data.selOpts['ADM'] = {} #{'thres': 0.01, 'inds': {'Type':'L', 'Eke':1.1}}
# data.setSubset(dataKey = 'ADM', dataType = 'ADM', sliceParams = {'t':[4, 5, 4]})
# data.ADMplot(keys = 'subset')
Aligned frame angular distributions¶
Finally, the observable angular patterns can be investigated, by looking at $|^{AF}T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\hat{k}_{L})|^2$.
NOTE data types for plots:
- 'r': real
- 'i': imaginary
- 'a': abs values
- 'phaseUW': phase, unwrapped (by energy)
# Run AF numeric code with unity AKQS parameters
data.afpadNumeric(AKQS = AKQS2pulseSub)
# Plot panels with padPlot
Erange = [1,40,1]
# This currently doesn't support the native dims, so reset these first and stack to standard LM representation (includes sum over MF terms)
for key in data.data.keys():
data.data[key]['TXafSum'] = data.data[key]['TXaf'].sum(['m','mu']).rename({'Lambda':'m'}).stack({'BLM':['l','m']})
data.data[key]['TXafSum'].attrs = data.data[key]['TXaf'].attrs
data.padPlot(Etype = 'Eke', Erange = Erange, dataType='TXafSum', pType='a2',
backend='mpl', sumDims=['Sym'], facetDims=['t', 'Eke'], selDims={'mu0':0},
reducePhi='sum', pStyle='grid', returnFlag=True) #, squeeze=True)
# TO TEST: should be able to add col_wrap=3 for Xarray plotter?
# See https://docs.xarray.dev/en/latest/generated/xarray.DataArray.plot.html
Differential case¶
# Set data
for key in data.data.keys():
for dataType in ['TXaf','TXafUnaligned']:
sumData = dataType + 'Sum'
data.data[key][sumData] = data.data[key][dataType].sum(['m','mu']).rename({'Lambda':'m'}).stack({'BLM':['l','m']})
data.data[key][sumData].attrs = data.data[key][dataType].attrs
# Set distributions
for pType in ['r','i','a','phaseUW']:
if 't' in data.data[key][sumData].squeeze().dims:
facetDims=['t', 'Eke']
else:
facetDims=['Eke']
data.padPlot(Etype = 'Eke', Erange = Erange, dataType=sumData, pType=pType,
backend='mpl', sumDims=['Sym'], facetDims=facetDims, selDims={'mu0':0},
reducePhi='sum', pStyle='grid', returnFlag=True, plotFlag=False); #, squeeze=True)
# Set plots by type (currently NOT done above!)
# if not data.data[key]['plots'][sumData]:
# data.data[key]['plots'][sumData][pType] = {}
# print(f"Setting data.data[{key}]['plots'][{sumData}][{pType}]")
# data.data[key]['plots'][sumData][pType]['pData'] = data.data[key]['plots'][sumData]['pData'].copy()
# data.data[key]['plots'][sumData][pType]['grid'] = data.data[key]['plots'][sumData]['grid'].copy()
# Dataout
outData = sumData + 'Out'
if not outData in data.data[key].keys():
data.data[key][outData] = {}
print(f"Setting data.data[{key}][{outData}][{pType}]")
data.data[key][outData][pType] = {'pData': data.data[key]['plots'][sumData]['pData'].copy()}
# Subtract
# data.data[key]['TXafUnalignedSum'] = data.data[key]['TXaf'].sum(['m','mu']).rename({'Lambda':'m'}).stack({'BLM':['l','m']})
# # Unaligned ref:
data.padPlot(Etype = 'Eke', Erange = Erange, dataType=sumData, pType='a',
backend='mpl', sumDims=['Sym'], facetDims=['Eke'], selDims={'mu0':0},
reducePhi='sum', pStyle='grid', returnFlag=True, plotFlag=False);
data.padPlot(Etype = 'Eke', Erange = Erange, dataType=sumData, pType='phaseUW',
backend='mpl', sumDims=['Sym'], facetDims=['Eke'], selDims={'mu0':0},
reducePhi='sum', pStyle='grid', returnFlag=True, plotFlag=False);
# Subtract distributions
# Note this will also define the sign of the phases
subData = 'sub'
data.data[key][subData] = {}
for pType in ['r','i','a','phaseUW']:
data.data[key][subData][pType] = {}
data.data[key]['sub'][pType]['pData'] = data.data[key]['TXafUnalignedSumOut'][pType]['pData'] - \
data.data[key]['TXafSumOut'][pType]['pData']
(E,Theta) differential plots¶
These should be respresentative of measurements...
NOTE data types:
- 'r': real
- 'i': imaginary
- 'a': abs values
- 'phaseUW': phase, unwrapped (by energy)
# PLOTS
pType = 'r'
data.data[key][subData][pType]['pData'].real.plot(x='Theta',y='Eke',col='t', col_wrap=3)
# PLOTS
pType = 'i'
data.data[key][subData][pType]['pData'].real.plot(x='Theta',y='Eke',col='t', col_wrap=3)
# PLOTS
pType = 'a'
data.data[key][subData][pType]['pData'].real.plot(x='Theta',y='Eke',col='t', col_wrap=3)
# PLOTS
pType = 'phaseUW'
data.data[key][subData][pType]['pData'].real.plot(x='Theta',y='Eke',col='t', col_wrap=3)
E-only diff plots - Eke scale¶
Summed over angular coord...
TODO: should pass full complex form here before sum!
Current case maybe incorrect - 10^-13 for r/i components? (Two anti-phase theta components cancel!)
# Sum over theta...
pType = 'r'
data.data[key][subData][pType]['pData'].sum('Theta').real.plot(y='Eke') #,col='t', col_wrap=3)
plt.title(f"Data type: {pType}")
# Sum over theta...
pType = 'i'
data.data[key][subData][pType]['pData'].sum('Theta').real.plot(y='Eke') #,col='t', col_wrap=3)
plt.title(f"Data type: {pType}")
# Sum over theta...
pType = 'a'
data.data[key][subData][pType]['pData'].sum('Theta').real.plot(y='Eke') #,col='t', col_wrap=3)
plt.title(f"Data type: {pType}")
# Sum over theta...
pType = 'phaseUW'
data.data[key][subData][pType]['pData'].sum('Theta').real.plot(y='Eke') #,col='t', col_wrap=3)
plt.title(f"Data type: {pType}")
THINK this should roughly compare to the results in Fig. 8.7 of Jan Tross' thesis - need to convert to harmonics scale, but he sees large feature around 4 - 5ps, flips sign between H9 and H13.
E-only diff plots - Harmonics scale¶
Summed over angular coord...
TODO: should pass full complex form here before sum!
Current case maybe incorrect - 10^-13 for r/i components? (Two anti-phase theta components cancel!)
data.data['orb5']['XS'].Ehv
# TODO: util function for this!
# See ep.util.conversion.conv_ev_nm
# # Define constants from scipy.constants
# h = scipy.constants.h
# c = scipy.constants.c
# evJ = scipy.constants.physical_constants['electron volt-joule relationship'][0]
# # Define output units - wavelength in m
# waveConv = 1e-9
# dataOut = (h * c)/(data * evJ)/waveConv
import scipy.constants
def setHHG(dataIn, w=785, Etype='hv', Eshift=0, dp=2):
"""
Set harmonics based on driving field photon energy.
Assumes Xarray input, and sets according to dataIn.Ehv/photonE.
Set a shift if required (since dataIn.Ehv may be off from real values).
dp is used to round values in output.
"""
# h = scipy.constants.h
# c = scipy.constants.c
# evJ = scipy.constants.physical_constants['electron volt-joule relationship'][0]
# waveConv = 1e-9
# # Photon energy
# photonE = (h * c)/(w * evJ)/waveConv
photonE = ep.util.conversion.conv_ev_nm(w)
print(f"Converting to harmonic E scale with \Omega={w} nm, photon E={np.round(photonE,3)} eV from Ehv, Eshift={Eshift}")
# Convert data scale to harmonics of specified E
dataIn['HH'] = ((dataIn.Ehv+Eshift)/photonE).round(dp)
dataIn.attrs['HHGconv'] = {'HH':dataIn['HH'].data.copy(),
'photonE':photonE,
'\omega':w}
print(f"Photon Erange [{dataIn.Ehv.data[0]}:{dataIn.Ehv.data[-1]}] set to harmonics \
[{np.round(dataIn.HH.data[0],2)}:{np.round(dataIn.HH.data[-1],2)}]")
# return testCoords
# testCoords = data.data['orb5']['XS'].Ehv.copy()
# setHHG(testCoords)
# data.data['orb5']['XS'].Ehv
# Set HHG axis in data.
# Includes Eshift to match 1st IP (14.5 eV), vs. 17.5 in orb5 results (but may be wrong orb!!!).
# TODO - check data attrs, FegeEng vs. orbIP
# setHHG(data.data['orb5']['XS'], Eshift=-3)
for pType in data.data[key][subData].keys():
setHHG(data.data[key][subData][pType]['pData'], Eshift=-4.2)
# Sum over theta... REPLOT ON HH SCALE
pType = 'phaseUW'
# pType = 'a'
# data.data[key][subData][pType]['pData'].swap_dims({'Eke':'HH'}).sum('Theta').real.plot(y='HH')
data.data[key][subData][pType]['pData'].sum('Theta').real.plot(y='HH')
plt.title(f"Data type: {pType}")
# Pull discrete harmonics...
HHaxis = np.arange(9,21,2)
data.data[key][subData][pType]['pData'].swap_dims({'Eke':'HH'}).sel(HH=HHaxis, method="nearest") \
.sum('Theta').real.plot(x='t')
plt.title(f"Data type: {pType}")
Comparison with expt¶
# Line plot
data.data[key][subData][pType]['pData'].swap_dims({'Eke':'HH'}).sel(HH=HHaxis, method="nearest") \
.sum('Theta').real.plot.line(x='t')
plt.title(f"Data type: {pType}")
Compare phases to Jan's thesis, Fig. 8.7.
- Similar behaviour.
- Switch in sign of phases between H9 and H13.
- t-axis may be slightly different in these cases (calcs herein are for 2-pulse alignment case).
Looks reasonable, given that the calcs herein use slightly different alignment parameters, and are for the HOMO only.
from IPython.display import Image
from IPython.core.display import HTML
# Image(url= "N2_phases_tross_thesis_fig8.7_crop.png", width=500, height=100)
Image(url= "N2_phases_tross_thesis_fig8.7.png", width=800, height=100)
Formalism¶
(TO CHECK - conj D terms...?)
The usual ePS eqns. for the MF results express the MF wavefunction as an expansion in the radial matrix elements plus some angular terms:
\begin{equation} T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})=\sum_{l,m,\mu}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)Y_{lm}^{*}(\theta_{\hat{k}},\phi_{\hat{k}})D_{\mu,\mu_{0}}^{1}(R_{\hat{n}}) \end{equation}The MFPAD is then determined from the square of this function. Note here that the polarisation term is usually expressed in the LF, then rotated into the MF.
The MFPAD can also be expressed in terms of $\beta_{LM}$ parameters, derived directly from the analytical square of the wavefunction. In this case the various angular momenta are coupled directly, and the observable is expressed as a set of $\beta_{LM}$s. This is usually preferable for computation of observables, however, since these are derived from the coherent square of the wavefunction, there is no phase information, so for coupling into further calculations one may prefer to work with the effective matrix elements or wavefunctions directly.
For the LF the same considerations apply, although one usually works with the $\beta_{LM}$s. For the LF the wavefunctions can be written as:
\begin{equation} ^{LF}T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\hat{k}_{L},R_{\hat{n}})=\sum_{l,m,\mu,\Lambda}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)D_{m\Lambda}^{l*}(\Omega)Y_{l\Lambda}^{*}(\hat{k}_{L})D_{\mu,\mu_{0}}^{1*}(\Omega) \end{equation}where the notation has been simplified (following [1]), and an additional term $D_{m\Lambda}^{l}(\Omega)$ (where $\Omega$ denotes Euler angles $(\phi,\theta,\chi)$) rotates the MF projection terms $m$ into the LF, and all Euler angles are integrated over (i.e. all molecular alignments/MF polarization orientations). Note this is usually written with the opposite labels, i.e. $m$ is the LF projection term and $\lambda$ the MF term, but we'll keep the ePS notation for consistency here.
For the AF, the distribution is essentially a convolution (assuming that the rotational wavefunction is separable). (See the AFBLM pages for more.)
ROUGHLY:
\begin{equation} ^{AF}T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\hat{k}_{L},R_{\hat{n}})=\intop d\Omega\sum_{K,Q,S}\sum_{l,m,\mu,\Lambda}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)D_{m\Lambda}^{l*}(\Omega)Y_{l\Lambda}^{*}(\hat{k}_{M})D_{\mu,\mu_{0}}^{1*}(\Omega)A_{Q,S}^{K}D_{Q,S}^{K*}(\Omega) \end{equation}Where $A_{Q,S}^{K}D_{Q,S}^{K}(\Omega)$ describes the axis distribution function in the LF Euler angles.
Here there is a triple product in $D$ - convert to non-conj terms (unnecessary?) then apply eqn. 3.118 in Zare:
\begin{eqnarray} \intop d\Omega D_{m\Lambda}^{l*}(\Omega)D_{\mu,\mu_{0}}^{1*}(\Omega)D_{Q,S}^{K*}(\Omega) & = & \intop d\Omega(-1)^{m-\Lambda}(-1)^{\mu-\mu_{0}}(-1)^{Q-S}D_{-m,-\Lambda}^{l}(\Omega)D_{-\mu,-\mu_{0}}^{1}(\Omega)D_{-Q,-S}^{K}(\Omega)\\ & = & 8\pi^{2}(-1)^{m-\Lambda}(-1)^{\mu-\mu_{0}}(-1)^{Q-S}\left(\begin{array}{ccc} l & 1 & K\\ -m & -\mu & -Q \end{array}\right)\left(\begin{array}{ccc} l & 1 & K\\ -\Lambda & -\mu_{0} & -S \end{array}\right) \end{eqnarray}Hence:
\begin{equation} ^{AF}T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\hat{k}_{L})=8\pi^{2}\sum_{K,Q,S}\sum_{l,m,\mu,\Lambda}A_{Q,S}^{K}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)(-1)^{m-\Lambda}(-1)^{\mu-\mu_{0}}(-1)^{Q-S}\left(\begin{array}{ccc} l & 1 & K\\ -m & -\mu & -Q \end{array}\right)\left(\begin{array}{ccc} l & 1 & K\\ -\Lambda & -\mu_{0} & -S \end{array}\right)Y_{l\Lambda}^{*}(\hat{k}_{M}) \end{equation}For the geometric computations, define this similar to existing AF $\Delta_{L,M}(K,Q,S)$ term (for AF $\beta_{LM}$ calculations):
\begin{equation} ^{AF}\Delta_{l,m}(K,Q,S)=(-1)^{m-\Lambda}(-1)^{\mu-\mu_{0}}(-1)^{Q-S}\left(\begin{array}{ccc} l & 1 & K\\ -m & -\mu & -Q \end{array}\right)\left(\begin{array}{ccc} l & 1 & K\\ -\Lambda & -\mu_{0} & -S \end{array}\right) \end{equation}Hence:
\begin{equation} ^{AF}T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\hat{k}_{L})=8\pi^{2}\sum_{K,Q,S}\sum_{l,m,\mu,\Lambda}A_{Q,S}^{K}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E){}^{AF}\Delta_{l,m}(K,Q,S)Y_{l\Lambda}^{*}(\hat{k}_{M}) \end{equation}Where $^{AF}\Delta_{l,m}(K,Q,S)$ is a channel coupling parameter.
[1] Underwood, Jonathan G., and Katharine L. Reid. “Time-Resolved Photoelectron Angular Distributions as a Probe of Intramolecular Dynamics: Connecting the Molecular Frame and the Laboratory Frame.” The Journal of Chemical Physics 113, no. 3 (2000): 1067. https://doi.org/10.1063/1.481918.
NEXT:
- Check result above more carefully, may have messed up some phases (likely). Should be non-conj term in $\mu,\mu_0$?
- Test numerically - should be able to compare numerical vs. beta results, similarly to LF test cases
Implementation - working from http://localhost:8888/lab/tree/dev/ePSproc/geometric_method_dev_2020/geometric_method_dev_pt3_AFBLM_090620_010920_161020_dev.ipynb to develop suitable geometric functions.
Versions¶
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter'])
# Check current Git commit for local ePSproc version
from pathlib import Path
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1
# Check current remote commits
!git ls-remote --heads https://github.com/phockett/ePSproc